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I. Introduction 


In this report, we summarize work done under NASA Contract 
No. NASl-15894 on advanced stability theory analyses for laminar 
flow control. The report consists of six sections, of which the 
last five are independent reports that summarize our progress on 
different aspects of this work. 

In Section II, we present a summary of our work on the 
SALLY stability analysis code for compressible flow problems. 

The resulting computer code seems to be at least an order of 
magnitude more efficient than previously developed codes. 

In Section III, we present a comparison of methods for 
prediction of transition using the incompressible SALLY computer 
code. It is shown that transition prediction by the envelope 
method and a new modified wave packet method are comparable 
in reliability but that the envelope method is more efficient 
computationally . 

In Section IV, we present a study of instability and. 
transition in rotating disk flow in which the effects of Coriolis 
forces and streamwise curvature on transition are investigated. 
Good agreement between the theory and experiments performed at 
NASA Langley Research Center has been achieved using the e^ method 
with N of order 10 at the onset of transition. 
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In Section V, we present a new 'linear' three-dimensional 
instability mechanism that predicts Reynolds numbers for 
transition to turbulence in planar shear flows that agrees well 
with experiment. We have extended the SALLY stability codes to 
compute directly this new instability mechanism. The results are 
in good agreement v/ith experiment. Subcritical transitional 
Reynolds numbers of order 1000 in plane Poiseuille and plane 
Couette flow and of order 2000 in pipe Poiseuille flow have 
been found. 

In Section VI, we present results obtianed using our 
stability analysis codes to study the finite-amplitude (nonlinear 
stability of axisymmetric pipe Poiseuille flow. Our results 
are in disagreement v/ith the earlier analytical work of Davey, 
Itoh, and Stuart. There appear to be no finite-amplitude 
nonlinear axisymmetric instabilities, in contrast to the 
predictions of the above-mentioned authors. The conclusion 
is that extreme care must be exercised in the application 
of Landau-Stuart-Watson perturbation ideas to such flows. 
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II . Efficient CoTnputi^t.i'ion o£ CoTnpressible Flow Stability of 
Three-Dimensional Boundary Layers 


In this section we present a method for calculati.on of com- 
pressible flow stability of three-dimensional boundary layers. 

The method is based upon a two point boimdary value (direct) 
approach and is more efficient than the commonly used shooting 
methods (1, 2) by an order of magnitude. 

The computer code SALLY (3) which was developed for LFC 
(laminar flow control) design applications employed incom- 
pressible stability theory. It is highly desirable to develop 
a COMPRESSIBLE SALLY code which can be used as a design tool 
for LFC applications. One basic prerequisite for such a 
code is a fast eigenvalue solver for compressible stability 
equations which comprise of a system of eight first order 
equations. The existing numerical methods generally use 
Runge-Kutta integration procedure along with Grara-Schmidt 
orthonormalization to control the parasitic error growth. 

This procedure is often very slow and thus very inefficient 
if used in a black box stability code such as SALLY. We de- 
scribe here an efficient method for obtaining eigenvalues 
of the compressible stability equations. This method will be 
used in the compressible version of SALLY. 

We write the governing equations as 

(A + B D + C) 4) = 0 (1) 

where cj) is a 5 row vector defined as 



( 2 ) 
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and A, B, C are 5X5 matrices with the following structure: 



In the above D = , where y is the normal boundary layeor 

coordinate and v' is the corresponding perturbation amplitude, 
u^ and w^ are the perturbation amplitudes in x and z directions 
and a, 3 are the respective wave numbers. T and P are the 
amplitudes of temperature and pressure fluctuations respectively. 

The boundary conditions for Equation (1) are 


4 


y = 0; 


0 


( 4 ) 


'f’l = ‘*>2 = S " '*’4 " 

y -V «; <t>2^ ‘t’3, *4 ° 

In order to solve eigenvalue problems posed by Equations 
(1) and (4) , we represent Equation (1) by central differences 
which results in a block-tridiagonal system of equations with 
5X5 blocks. The system of equations is solved using Lu 
factorization. We use inverse iteration procedure (4) for 
eigenvalue search. This method is very effective once a 
good guess for the eigenvalue is available because the con- 
vergence is cubic. 

If we write Equation (1) as 

L (|) = 0 

then, the inverse iteration algorithm for obtaining eigenvalue 
can be written as 

(L - \ I) ‘t’K + 1 = f -t'k 

(L - X,^ \ + 1 = fk 

^k + l~ '*'k + l^‘*’k + l y'*'k + l‘*’k + l 

The finite-difference method presented above is second 
order accurate. However, the accuracy of the eigenvalue ob- 
tained by this procedure can be increased by Richardson's 
extrapolation. For this purpose we use Neville’s algorithm. 

If 

= X (h^) (i = 0, m) 
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j = 1, . . . m 

X “ 0 / • • • m““ ^ 

(6) 

We present here eigenvalue results for the leading edge 

region (R = 150) of a 35° swept back wing of infinite span. 

Table 1 gives the eigenvalues calculated for different grids 

= — ) and the extrapolated values. The sequence for 

^ i 

h chosen is that proposed by Bulirsch and Stoer (5) which 
i 

gives an eigenvalue converged to five significant digits. 

The cost of computation obviously depends upon the required 
accuracy. Generally, an eigenvalue accurate to 3 significant 
digits can be obtained in less than 2 seconds on Cyber 175. 
Figure 1 shows eigen-function distribution for the same flow 

as above. 

Another way to increase accuracy is to use Chebyshev 
spectral method (6). The spectral methods provide infinite 
order of accuracy and are extremely desirable particularly 
for high Reynolds number applications. Let us represent 

Equation (1) as 

L ; = f 

sp 

where L is a spectral operator, 
sp 

The direct solution of (7) by Gauss elimination methods 

3 

would require order (N = number of points) and order N 


where h^ is the grid size, then 


p(j) _ + 

1 + 1 


1-1 


p; 

1 


j) 


hf \ 2 


- 1 
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arithmatic operations. We describe here a method which permits 
solution of (7) using order N storage locations with the number 
of arithmetic operations of the order of N log N. Specifically 
we use Chebyshev acceleration scheme (7) : 

_ .n + 1 _ ,(n) . /I \ A \ 
ap ^ ap ) n ^ ' n ^ ^ 


- a <. 1 ^ (LspU'"^ - f) (8) 

where is an approximate finite-difference operator and 

is a relaxation parameter. The error in the solution of 
Equation (7) is decreased by a factor of 10^ after about 
9 iterations of Equation (8). 


0) 

n 
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TABLE 1 






RICHARDSON 

EXTRAPOLATED 

EIGENVALUES 

X— (a = .2724 

, 3 = “ .2915) 



i 


m = 


1 


2 


3 

4 


5 


20 

1 

-.2541052E-01 

.662266E-02 

-.2505890E-01 

.6155269E-02 






! 

40 

1 

*.2514681E~01 

.6272069E-02 

-.2508481E-01 

.6158182E-02 

-.2508804E-01 

615854GE-02 




. 1 

60 

9 

-.2511236E-01 

.6208798E-02 



-.2508700E-01 

•6158786E-02 

- . 2500698E-01 .6158802E-02 





■ 



-.2508645E-01 

.6153633E-02 



-.2508797E-01 .6159012E-02 

-.2508800E-01 

.6159018E-O2 


80 

1 

-.2510183E-01 

.6186852E-02 



-.2508786E-01 

.6158987E-02 




25087 85E -01 .6159061E-02 

120 

5 

-.2509352E-01 

.6171322E-02 

-.2508751E-01 

.6153899E-02 

-.2508786E-01 

.6150047E-02 

-.2508786E-01 .6159057E-02 

.2508786E-01 

.6159 060E-02 

j 

I 

j 





-.2503777E-01 

.6159010E-02 






1 

160 

6 

-.2509100E-01 

.6165936E-02 










VD 
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/ 
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III. Comparison of Methods for Prediction of Transition by 
Stability Analysis 

ABSTOACT 

Several methods of transition prediction by 
linear stability analysis are compared. The spectral 
stability analysis code SALLY is used to analyze 
flows over laminar flow control wings. It is shown 
that transition prediction by the envelope method 
and a new modified wave packet method are comparable 
in reliability but that the envelope method is more 
efficient computationally. 
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NOMENCLATURE 


A 


c 

f 

k 

L 

N 

R 

Re 

c 

s 


n 


t 

U 

u. 



w 

w 

X 

y 


z 


insxiitiiim disturbance amplitude 
Chebyshev coefficients 
wing chord 

dimensional frequency 
wave vector 

algebraic mapping parameter 
N-factor - £n A/A^ 

* 

displacement thickness Reynolds number, U 6 /v 

chord Reynolds number, c/’^oo 

arc length along an arbitrary path on the wing 

Chebyshev polynomial 

time 

unperturbed x-velocity in the boundary layer 
potential flow vector at edge of boundary layer 
x-component of 

incoming free stream velocity 
group velocity vector 

perturbation velocity in the y-disection 
unperturbed z-velocity in the boundary layer 
mapped coordinate normal to wing surface 
coordinate in the direction of the normal chord 
coordinate normal to the wing surface 
coordinate along the wing span 
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x-wave number 
angle of attack 
z-wave niunber 
frequency 

displacement thickness 
wave length 
kinematic viscosity 
wing sweep angle 

angle formed by the wave niimber vector with the 
x-axis 

angle formed by the group velocity vector with 
the x-axis 

angle formed by the potential flow vector with 
the x-axis 

eigenfunction; defined in Eq. (3) . 



INTRODUCTION 


In this section, several methods of transition prediction 

using linear stability analysis are compared. The incompressible 
linear stability computer code SALLY is used in various ways 
to study three-dimensional boundary layer flow over laminar 
flow control (LFC) wings. Here we compare the so called 
envelope method with wave-packet methods , to predict 
transition. We conclude that the envelope method is at 
least as reliable as the more complicated and less efficient 
wave packet method. 

Consider the stability of three dimensional laminar 
flow over swept wings with sweep angle 0 . The coordinate 
system used on the wing is depicted in Fig. 1. The x-axis 
is in the direction of the noimial chord, the y-axis is 
normal to the surface of the wing while the z-axis is 
along its* span. 

Neglecting the curvature of the wing surface, 
compressibility effects, and non-parallel flow effects, 
linear disturbances satisfy the Orr-Sommerfeld equation 


, d^ 2 .2,-2 . 

( — 9 - a - 3 ) <|5 

dy^ 


2 2 2 

* iR{ (aU + 3W - 03 ) [-^ - - 3^]c|) - (D 


dy 


dy 


dy 
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with the boundary conditions 

4>(0) = ^ (0) = 0; 4) (oo) bounded. (2) 

Here the perturbation velocity in the y-direction is assumed 
to be of the form 

V = , (3) 

U(y) and W(y) are the (unperturbed) laminar boundary 
layer velocities in the x- and z~directions, respectively, 
and R is the Reynolds number. It is assumed that all 
variables are non-dimensionalized with boundary layer 
scaling. 

Equations (1) - (3) constitute an eigenvalue 
problem for the frequency o) and wavenumbers a, 3 - 
For given Reynolds number R, this eigenvalue problem 
provides a complex dispersion relation of the form 

0) = (o(a,3) (4) 

relating the complex parameters a, 3 and o). 
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Semi-empirical methods to predict transition on 

LFC wings are based on tracing the evolution of modes 

1 

across the wing. An appropriate N- factor for transition 
correlation is defined as the (logarithm of the) total 
growth factor across the wing (see below) . A good transition 
predictor is one for which transition occurs at nearly constant 
N for a wide variety of wings and flow conditions. 

For natural transition, disturbances of all frequencies 
are present on the wing surface. In this case, there are 
many optional ways to compute N factors. The first choice 
is between temporal and spatial stability theory. In 
temporal theoiry, a and B are real while o) is complex? 
the mode grows in time if Im(oj) >0, but the mode does 
not grow in space. An N-f actor for transition correlation 
may be defined as 

- g 

N = I Im(o))/|Re (v ) I ds 

where v^ =Oo3/3a, 3o)/3B) is the (complex) group velocity 
and s is the arclength along an appropriate curve on the wing. 
The N-f actor (5) is not fully defined until a prescription 
is given for singling out a specific mode at each position ^ 
on the wing and for defining a specific curve on which, to 
integrate. We shall return to these questions in Sec. 2. 

In spatial stability theory, w is real but a and/or 
B may be complex. Again, there is arbitrariness in the 
definition of an appropriate N-factor because of the variety 
of excitable modes on the wing. 
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WAVE PROPAGATION IN BOUNDARY LAYERS 


The complex eigenvalue relation (4) provides 
two. real relations among the three complex quantities 
tt/3, and co. In temporal stability theory, the 
requirements that a and 3 be real provide two more 
conditions so there remain two arbitrary parameters 
among Re(a), Re(3)/ Re(co), and Im(o)) . 

There are several ways to remove this arbitrainess 
in the computation of the growth factors N. In the 
envelope method^, Im(u>) is maximized with respect 
to a at fixed Re o) [which then determines a, 3 
and 0 ) uniquely at each point on the wing] and the 
curve in (5) is defined to be everywhere tangent 
to Re (V ) . 

y 

With spatial stability theory, there remains 

/ ■■ ^ 

three independent real parameters among a , 3 

Re(w) once the eigenvalue condition (4) is satisfied. 

One possibility is to require that the direction of 

most rapid growth, which is parallel to the vector 

(-lm(a), -Im(3))/ be parallel to Re(v^) and that 

the resulting value of the most rapid growth rate 

be maximized with respect to the remaining two 

3 

independent parameters . 

Alternatively, it is possible to use wave 
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packet theory to remove the arbitrariness in the definition 

of N-f actors. For a conservative dynamical system, 

kinematic wave theory implies that a wave packet 

propagates in physical and wavevector space according 

4 

to the Hamilton-Jacobi equations 


dx 

— 

dt 

Soi 

dz 


dt 

33 

da 

_ 3m 

dt 

“ 3x 


= _ 

dt 

3z 


( 6 ) 

(7) 

(8) 

(9) 


2 

Nayfeh considered the extension of Eqs. (6)- 

(9)* to non-conservative systems where a# 3/ and m can 

be complex. Then, and may also be complex. 

For a physical solution with real x,z, and t to 

exist, (6) and (7) imply that the group velocity 

(oj ,Wo) must be real. Nayfeh proposed the computation 
a p ^ 

of wave packet solutions determined by the six independent 
conditions: (i) the eigenvalue condition (4) ; 

(ii) Im 0 ) = Im 0 )^= 0; (iii) Re m fixed; (iv) 

a p , “ ■ 

Re 3 fixed; and (v) dx/dt = m , dy/dt = cjo* Under 

oc p 

these conditions the N-f actor is determined by 
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( 10 ) 


t 

N = / [-0) Im(a) •~a)^Im(3)+Im(a)) ]dt 

4. ^ p 

Finally we study a modified non-conservative wave 
packet formulation in which a,^, and oj are determined 
by; (i) the eigenvalue condition (4) ; (ii) 

Im = Im 0 )^ = 0; (iii) Re o) fixed with Im w = 0; and 

(iv) dx/dt - 0 )^, dy/dt = o)- . The motivation for these 

a p 

latter conditions is simply that laminar flow over a LFC 
wing may be assiamed steady so a wave packet should propagate 
at fixed real frequency. However, there is less justification 
for assuming Re 3 is fixed as in Nayfeh’s formulation, 
because the flow is not homogeneous in space. The N- 
f actor is given by (10) with Im(o)) = 0. 

With Nayfeh*s formulation of the wave packet equations, 
the growth factor N is a function of the independent 
variables Re u) and Re 3 , while N is a function of only 
Re 0 ) in our wave packet formulation. Therefore, 
maximization of N over all allowable packets is computation- 
ally more efficient with our formulation. We have performed 
computations (not reported in detail here) with Nayfeh's 
wave packet formulation and have found the computations to 
be extremely sensitive, with realistic solutions satisfying 
the required constraints at some wing locations but not 
at others and the overall N factor at transition highly 
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valuable. We do not believe these latter effects 

originate in the numerical scheme. In any case, the conclusion 

of the present paper that the envelope method is at least 

as reliable as the wave packet method and that it is considerabj.y 


more efficient would not be changed by comparisons with 


results obtained by Nayfeh^s wave packet formulation. 


NUMERICAL METHOD 

In the computer code SALLY^, Eqs. (1) - (2) 

are solved using a spectral method based on Chebyshev 
. 5 

polynomials . The boundary layer direction y 
0 £ y < oo, is mapped into the finite interval 
-1 £ w < 1 by the algebraic mapping 


2 

y+L 


1 


( 11 ) 


and 4> (y) is approximated as the finite Chebyshev 
polynomial series 

= 1 T„(w) (12) 

n=0 " “ 

The resulting algebraic eigenvalue problem is solved 
globally (if a guess for the eigenvalue is not 
available) by a generalized QR algorithm or locally 
(if a good guess is available) by inverse Rayleigh 

g 

iteration . The resulting scheme is very efficient 
and accurate. 
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The properties of the laminar boundary layer 
profiles required to solve (1) - (2) are obtained using 
a compressible boundary layer code for swept tapered 
wings developed by Kaups and Cebeci^’ 

The code SALLY also performs a number of optional 
computations, including: (i) computation of maximum 

amplification among all wavelengths and propagation 
angles; (ii) computation of amplification at fixed 
frequency and fixed wave length; (iii) computation 
of amplification at fixed frequency and fixed propagation 
angle; (iv) computation of maximum amplification at 
fixed frequency; and (v) computation of wave packet 
solutions according to the prescriptions discussed 



RESULTS 


7 

Burrows has reported flight transition data 
taken at Cranfield for a large, untapered, 45® swept 
half wing mounted as a dorsal fin upon the mid-upper 
fuselage of an Avro Lancaster airplane. The airfoil 
section was made-up of two semi-ellipses, one of which 
constituted a faired trailing edge and the other 
corresponding to the leading edge of a 10 percent 
thick airfoil, with effective chord of 10.83 feet, 
measured in the free stream direction. The location 
of the beginning of transition in the Cranfield data 
was estimated as given in Ref. 8 . Two of the Cranfield 
flight rests were chosen for correlating transition 
using wave packet theory. 

In the first test case, calculations were made 

g 

for a chord Reynolds number of 11.7 x 10 and -2® 
angle of attack. In this flow, transition begins at 
x/c = 5.5%. A maximum N factor of 7.6 was obtained 
at a frequency of 1250 Hz both with the envelope 
method and the modified wave packet method. 

The predicted variation of the N factor up 
to the transition location was almost identical for the 
envelope method and the modified wave packet method. 

We also compute the solution of the conservative wave 
packet equations (6) - (9) in which only the real 
parts of equations (6) - (7) are taken while (8) -(9) 
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are solved in their full complex form. The resulting 
N factor at transition is 5-2. The variation of N 
factor with x/c for the various methods is plotted 
in Fig. 2. 

Wave angle , wave length and the direction of 
the group velocity as predicted by the envelope and 
wave packet methods are given in Figs. 3-5. Although 
the results are qualitatively similar, there is 
appreciable quantitative difference in these parameters 
at the transition location. It is surprising that the 
N factor calculated by the envelope and modified wave 
packet methods are the same. 

In the second test case, the angle of attack 
of the wing was changed to zero. In this case, 
transition occurred experimentally at x/c = 7%. 

The envelope method gave an N factor of 10.8 at 
a frequency of 1000 Hz. The wave packet method gave 
a maximum N factor of 10.5 at a frequency of 
1200 Hz, which is close to the prediction of the 
envelope method. The variation of N factor with 
x/c is plotted in Fig. 6. The predictions of the 
conservative wave packet approximation and a fixed 
wavelength, fixed frequency integration are also 
plotted in this figure. The conservative wave packet 
approximation gave an N factor at transition of 
8.6 rather than 10.5. 



on N 


Figure 7 shows the influence of frequency 
factor at transition for the wing as 
predicted by the wave packet theory. Wave angle, 
wave length and direction of the group velocity 
for this particular wing are shown in Figs. 8-10. 
Again there is substantial quantitative difference in 
the predictions of the two methods. 


27 


CONCLUSIONS 

Calculations were made for a Cranfield 45° 

6 

swept wing with Re = 11.7 x 10 using a modified 
wave packet method and the envelope method. Both 
methods gave an N factor of 7.6 at transition 
location for an angle of attack, = - 2°. For 

= 0°, the envelope and modified wave packet 
methods gave N factors of 10.8 and 10.5, 
respectively. Since it may be argued that the 
wave packet method is physically more relevant for 
predicting transition in three dimensional boundary 
layers, it was initially hoped that the wave packet 
method might give more consistent transition N 
factors. However, the results show that the wave 
packet method provides N factors which are at best 
as consistent as those of envelope method. Since the 
wave packet method is at least 3 times as expensive 
to use as the envelope method, the latter is recommended 
for engineering design calculations. 
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FIGURE CAPTIONS 


Figure 1. A plot of the coordin£.te system on a swept 
wing . 

Figure 2. A plot of N versus percent of chord 
x/c for various methods applied to a swept wing 
of an Avro Lancaster airplane at -2® angle of 
attack. Solid curve; modified wave packet method 
and envelope method at f = 1250 Hz which gives nearly 
the maximum N at the transition point. Dashed 
curve; result of integrating equations (6) - (9) 
across the wing with (6) and (7) replaced by 
their real parts. The curves are plotted from the 
beginning of the unstable flow region until the 
transition point at x/c =5.5%. 

Figure 3. A plot of wave propagation angle versus 
x/c for the same flow as in Figure 2. 

Figure 4. A plot of wavelength versus x/c for the 
same flow as in Figure 2. 

Figure 5. A plot of the direction of the group velocity 
for the same flow as in Figure 2. 

Figure 6. Same as Fig. 2 except for the wing at 0® 
angle of attack. In addition to the results of the 
wave packet methods and envelope method, the N 
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factor obtained by integrating a fixed 
wavelength, fixed frequency mode across the wing is 
given. Here N is given by (5) and the mode is 
- detexrmined by the six real conditions: (i) Eq. (4) 

(ii) Im a = -Im 3 = 0; (iii) i/c- 0.001; (iv) 

Re 750 Hz. 

Figure 7. Variation of N at transition versus 
frequency obtained using the modified wave 
packet method for the same flow as in Figure 6. 

Figure 8. A plot of wave propagation angle versus 
x/c for the same flow as in Figure 6. 

Figure 9. A plot of wavelength versus x/c for the 
same flow as in Figure 6. 

Figure 10. A plot of the direction of the group 
velocity for the same flow as in Figure 6. 
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Figure 6. 
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Figure 7. 
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Figure 10. 
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Instability and Transition in Rotating Disk Flow 


ABSTRACT 


The stability of three-dimensional rotating disk 
flow is investigated, including the effects of Coriolis 
forces and streamline curvature. The results show that 
the critical Reynolds number for establishment of 

stationary vortex flow is Rq = -287.? These vortices 

spiral outwards at an angle of about 11.2° and 
transition to turbulence occurs when their total 
amplification is about e^?". We also report new 
experimental results on the spatial growth rates of 
the stationary vortices. It is shown that our 
analysis gives grov/th rates that compare much better 
with the experimental results than do results obtained 
using the Orr-Sommerfeld equation. Our calculations 
also indicate the existence of weakly unstable propagating 
(Type II) modes at low Reynolds numbers (R^ = 49) . 



1 . 


INTRODUCTION 


The prediction of transition in three dimensional 
boundary layers [1”3] is a subject of both fundamental 
and practical importance in fluid mechanics. Practical 
interest in the subject centers on the design of laminar 
flow control (LFC) wings that promise significant 
improvement in airplane fuel efficiency. At present, 
the most useful tool for transition prediction in such 
flows is the so-called e^ method [4]. Hefner and 
Bushnell [5] and Malik and Orszag [6] show that 
the exponent N (called the N factor) is of the 
order 7-11 when transition occurs on LFC swept 
wings. 

The instability mechanism exhibited in the 
leading edge region of a swept wing is similar to that 
found in the boundary layer on a rotating disk, since 
both have mean cross flow profiles with inflection 
points. More details on the similarities between the 
two flows is given in Ref. [7]. However, the rotating 
disk flow is more convenient to study in view of von 
Karman's exact steady solution of the Navier-Stokes 
equations [8] . 

Using hot-wire techniques , Smith [9] observed that 
sinusoidal disturbances appear in a rotating disk boundary 
layer at sufficiently large Reynolds number. About 32 
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oscillations were observed within a disk rotation period 
and analysis indicated that the disturbances propagate at 
an angle of about 14® relative to the outward drawn radius 
(where the direction of disk rotation defines positive 
angles). Later, in a remarkable study using the china- 
clay technique, Gregory et al [10] observed about 30 
vortices over the disk spiraling outwards at an angle 
of about 14'^ (that is, their normals make an angle of about 
14° with the outward dravm radius).. These vortices, which 
appeared stationary relative to the rotating frame of 
the disk, were first observed at a Reynolds number -430 

[where R is defined after (2.14) below]. Transition 
to turbulence occurred at The overall conclusions 

drawn by Gregory et al have been confirmed in later 
investigations [11,12]. It has been found that there 
are about 30 stationary vortices whose normals make an 
angle of 11° - 14° with the radius. However, there seems 
to be considerable controversy over the value of the 
critical Reynolds number which in our view can be attributed 
to the measurement techniques used. Further, at low 
Reynolds numbers, Fedorov et al [13] observed only 14-16. 
vortices with normals lying at an angle of about 20°. 

Stuart [10] analysed the linear, inviscid 
stability of rotating disk flow. However, the neglect of 
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viscosity resulted in the prediction of 113-140 vortices 
over the disk which is about four times larger than the 
observed value. Brown [14] extended Stuart's work to the 
viscous case by applying the Orr-Sommerf eld equation to disk flow. 

Using temporal instability theory. Brown found R sj 178, 

C 

which is much less than the observed value. Recently, 

Cebeci and Stewartson [3] solved the Orr-Sommerf eld 
equation for rotating disk orofiles using spatial stability 
theory and found - 170. They also suggested that 

wave packets propagate in three dimensional flows in such 
a way that dct/d3 is real. Using this condition, Cebeci 
and Stewartson correlated transition using the e^ method 
and found the N factor at transition 510) to be 

X 

about 20 which is much higher than that found for LFC 
wings ([5] , [6]) . Bushnell (private communication) argues that 
a higher N factor may be required for transition in disk flow than 
on LFC wings because the boundary layer is rotating with 
the disk while the external disturbances in the surroundings 
are not. Consequently, there is no direct coupling between 
the external disturbances and the instability waves in the ^ 
rotating disk boundary layer. 

The Orr-Sommerfeld stability equation neglects 
the effects of Coriolis forces, streamwise curvature, and 
nonparallel flow.' In Ekman layer flow, Lilly [15] has 
shown that the Coriolis force has a significant effect on 
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stability at low Reynolds numbers. Lilly showed that the 
critical Reynolds number for appearance of stationary 
vortices is higher (R^ - 115) when the Coriolis 
force is included in the analysis than when it is neglected 
(Rq ~ 85). In addition to the stationary vortices, 
he showed the existence of a parallel instability caused 
by the Coriolis force at much lower Reynolds numbers. 

Such an instability mechanism was also observed in the 
Ekman layer experiments of Faller and Kaylor [16] and 
Tatro and Mollo-Christensen [17] . The Ekman layer and 
the rotating disk are similar in that both are three- 
dimensional boundary layer flov/s in which rotation plays 
a significant role. Lilly's results suggest that the 
inclusion of the Coriolis force in the stability analysis 
of rotating disk flow may also lead to a higher critical 
Reynolds number for stationary vortices which is in better 
agreement with observations. 

In this section we present a stability analysis of 
rotating disk flow in which the effects of Coriolis force 
and streamline curvature are included. The resulting 
sixth order system is solved numerically by a Chebyshev 
spectral method [18-19]. We also follow the evolution 
of the disturbance modes using the envelope method 
[1,6] and calculate the N factor at transition. 


46 



The work of Kobayashi et al [12] , which appeared during 
the final stages of the present study, also includes the 
effects of the Coriolis force and streamline curvature. 

We will comment on this v/ork in Sec. 6. 

We also report results of an experimental investigation 
of rotating disk flow in which the growth rates of the 
boundary layer disturbances were measured using a hot-wire 
probe. While only limited experimental results are now 
available, the agreement with the present theoretical 
analysis seems encouraging. We know of no previous 
experimental study in which growth rates were measured. 
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2, THEORETICAL ANALYSIS 

Consider an infinite plane disk rotating about 
its axis with angular velocity Q . We take cylindrical 
coordinates r*,0,z* with z* - 0 being the plane of 
the disk. Assuming the fluid to lie in the half-space 
z* > 0, the governing dynamical equations are, in the 
rotating coordinate system, 

^ , V* 3u* , 9u* v*^ _ 1 9p* 

9t* 9r* r* 90 ^ 9z* ^ r* p dr* 


+ v( 
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+ u* — + — — + w* — - - - ^ 
3t* 3r* r* 39 '' 3z* p 32 * 
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2 2 
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— T 2 

30 3z* 


1 3w* . 
r 3r*^ 


(.2,3) 


3u* 3v* 3w* u * 

3r* r* 36 3z* r* 


(2,4) 


where (u*,v*,w*) are the velocity components in the 
(r*,6,z*) directions, respectively/ p* is the pressure, 
p is the density, and v is the kinematic viscosity. 

Von Karman ' s exact solution of C2.1) - (.2,4) for 
steady laminar rotating disk flow is obtained as follows 
[8] . Let u,v,w, and p denote the steady state values 
of u*,v*,w*, and p* respectively. Defining 

u = r*^^F(z) , V = r*fiG(z) , w -AjQH(z) , p =pvf^p(z) , (2.5) 

* 

where z=z /P/v then the Navier-Stokes equations (2.1) - 
(2.4) reduce to the following equations for F,G,H and 
P: 

F^ - (,G+1)^ + F'H - F" = 0 (2.6) 

2F(G+1) + G’ H - G" = 0 (2.7) 
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P' + HH 


H" 


0 


( 2 . 8 ) 


2F + H' = 0 


(2.9) 


where the prime denotes differentiation with respect to 
z. The boundary conditions are 


F=0, G-0, H=0 at z=0 


F - 0 , G = “ 1 for 2 ->■ a> 


( 2 . 10 ) 


Now we study the evolution of infinitesimal small 

disturbances imposed on the steady flow governed by 

Eqs . (2.5) - (2.9). Let r^ be the radial location near 

which the analysis is to be made. Using r Q as the 

e 

reference velocity, 6* = /v7fi as the reference length, 

2 2 

and pr^ Q as the reference pressure, the perturbed 
nondimensional velocities u,v,w and pressure p can 
be written as 


u(r,6,z,t) = ~ F(z) + u(r,0,z,t) 
v(r,e,z,t) = I G(z) + v(r,0,x,t) 
w(r,0,z,t) = ^ H(z) + w(r,0,z,t) 


( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 
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p(r,e ,z,t) 


(2.14) 


= -5- P (z) + p (r ,e ,z ,t) 

R 

Here the nondimensional radius is r == r*/TI7v, the Reynolds 
number is R = r and r^ corresponds to r - R, 

Substituting Eqs. (2.11) - (2.14) in the Navier- 
Stokes equations (2.1) ~ (2.4) and linearizing with 
respect to the perturbations gives : 
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(2.18) 


9ti , 1 3w ^ u 

3r r 30 3z r 

The boundary conditions are that U/V, and w vanish 
at z = 0 , 00 . 

For R>>1, the system (2.15) - (2.18) may be 

consistently approximated by replacing factors of r 

-2 

by R and neglecting terms of order R and smaller. 

The replacement of r by R at this stage of the 

calculation implies that we neglect some nonparallel flow 

effects. These effects are now under study. The neglect 

— 2 

of terms of order R and smaller has little effect on 
the results discussed below, as we verified by computations 
in which they were included. 

Replacing factors of r by R in (2.15) -- (2.18) 
gives a set of equations that are separable in r,0,t 
so that the perturbation quantities may be assumed to 
have the form 

(V^,^,w,p) = (f(z),h(z),(|,(2),Tr(z))e^'“^^®^®~“^’ (2.19) 

With this assumption, Eqs. (2.15) - (2.18) become (not yet 

— 2 

dropping termis of order R ) 

i(aF+aG-o))f+F> + iam = ^ [f’-X^f - Ff + 2(G+l)h - Hf] 

+ i [iaf - 2i3h]- f (2.20) 

r"^ r-^ 
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KaF + 3G-w)h + G'(f) + iBTr = i [h"~A^h~Fh-2 (G+1) f-Hh ’ ] 

+ [iah + 2i 3 f ] -- i h (2.21) 

R 

KaF + 3G-w) 4)+ tt’ ^ \ cc (p (2.22) 

R R/i 

(ia+|)f + ±3h + (|)' = 0 (2.23) 

where + 3^. 

Eliminating tt from (2.20) - (2.22) by means of 
(2.23) gives, neglecting terms of order and smaller, 

I [i(D^-A^) (D^-X^) + R(aF+3G-(D) (D^-X^) - R (Xf"+ 3G” ) -iHD (D^-X^) 

~ iH'(D^-X^) - iFD^] (|>+^[2.(G+l)D + 2G ' ] n = 0 (2.24) 

R 

[2(G+l)D-’iR(aG'-3F* ) ] (^ -+i[i{D^~X^) + R(aF+3G-w) - iHD - iF]n= 0 


(2.25) 

_ —2 — 2 

and where D = d/dz, a = a-i/R,X = aa+ 3 and 
n = ah “ 3f is proportional to the z-component of the 
perturbation vorticity. The final result (2.24) - (2.25) is a con- 
sistent set of stability equations valid to order R~^. 

The boundary conditions for the sixth order system 
(2.24) - (2.25) are 
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4^(0) = (|)' (0) = n(0) =0 

(2.26) 

(|) (oo) = (|) » (oo) = (oo) = 0 

Note that if the Coriolis force and streamline 
curvature effects are neglected, the above system 
reduces to the fourth-order Orr-Somir.erfeld equation; 

+ R{aF-(6G-u) - R (aF"+gG " ) ] <f) = 0 

(2.27) 

In Sec. 6, we report numerical results for both the 
sixth order system (2.24) - (2,25) and the fourth order 
equation (2 . 27) in order to study the effect of Coriolis 
force and streamline curvature terms on the stability of 
flow due to a rotating disk. 
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3 . NUMERICAL METHOD 


We solve the Orr-Sommerfeld equation (2.27) in 
the computer code SALLY [1] by using a spectral method 
based on Chebyshev polynomials [ 18-19 ] .Here we extend the 
method to solve the sixth order system (2.24) - (2.26). 

The boundary layer coordinate z, 0 _< z < co is 
mapped into the finite interval - 1 ^ ^ <1 by the 
algebraic mapping 




z+L 


1 


(3.1) 


where L is a scale parameter chosen to optimize the 
distribution of points in Then (z) and n (z) 

are approximated as the finite Chebyshev polynomial 
series 


M 


n=0 


(3.2) 


M 


n (z) = 


y b T (r) 
n=0 


(3.3) 


Substituting (3.2), (3.3) in (2.24) - (2.26) and 
collocating [19] at the discrete points = cos -iij/M 
(0 £ j £ M) gives the algebraic eigenvalue problem, 
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n 


n 


= 0)B 


n 


n 


(3.4) 
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where A and B are 2(Mfl) ^2(Mfl) matrices. The 
eigenvalue problem (3.4) is solved globally (if a 
guess for an eigenvalue is not available) by a generalized 
QR algorithm or locally (if a good guess is available) 
by inverse Rayleigh iteration [20]. The resulting scheme 
is very efficient and accurate. In the present 
calculations the optimum value of the scaling parameter 
L was found to be about 1,8. In most of the calculations 
reported below, M = 34 so 35 Chebyshev polynomials were 
used. 

The accuracy of the method was tested in several 
ways. First, the number of retained polynomials, M+ 1 
was varied to check the accuracy of the eigenvalues and 
eigenfunctions. Then, calculations were made for the 
stability of Ekman flow. Comparisons were made with the 
results obtained by Lilly [l5] for R = 65, 110, 

150, 300 and 500 with good agreement. 

For rotating disk flow, the global method gives only 
one unstable eigenvalue [Im(w) >0] for R > 150 that 
is insensitive to M. However, spurious unstable modes 
appear for lower R which are discarded as unphysical 
because they are very sensitive to M . 
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4.. TRANSITION PREDICTION USING THE METHOD 

In three dimensional flow, the dispersion relation 
is given by the complex relation 


0 ) = 0 ) (a , 6 ) 


(4.1) 


where and w are, in general, complex. Therefore, 

there are four arbitrary real parameters among a, 3 and 
0 ). There are several ways [6] to remove this arbitrariness 
In the present study we employ the envelope method [1] . 

Here the four conditions are obtained by using temporal 
stability theory [in which Im (a) = Im{g)- 0] and by 
maximizing Im(w) with respect to a, 3 at fixed Re(w). 

The N factor is then given by 


N = / 

s. 


Im(o3) 


I Re(Vg) I 


ds 


(4.2) 


where ~ (complex) group velocity and 

s is the arc length along the curve whose tangent is the 
real part of the group velocity. Noting that 


ds = 
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Eq. (4.2) can be written as 


N = 


/ 

R, 


T 


Im (oj ) 

Re ( CO ) 
a 


dR 


(4.3) 


Here the subscripts C and T indicate critical (linearly 
unstable) and onset of transition, respectively. 
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5. EXPERIMENTAL STUDY 

An experimental program was established to quantitatively 
study the flat rotating disk flow with particular attention 
given to measurement of the growth of boundary layer disturbances 
as a function of Reynolds number. The experimental set-up 
is described along with the measurement techniques. 

Measured mean velocity profiles are compared with the exact 
solution and an analysis of mechanical disk vibrations is 
presented. 

Rotating Disk Apparatus 

The experimental apparatus is shown schematically 
in Figure (1). A 457 mm diameter, 12.7 mm thick Plexi- 
glass disk was attached to a 50.8 mm diameter aluminum 
shaft by means of two parallel aluminum mounting disks. 

Shim stock inserted at various circumf rential locations between 
the mounting disks was used to control alignment of the 
Plexiglass disk and to compensate for asymetrical flexure 
of the Plexiglass introduced by mounting stresses. The 
drive shaft was inserted between two pre-mounted, self- 
aligning ball bearings and was driven by a 1/4 Hp, 1725 
RPM AC motor through a 2:1 belt and pulley speed reduction 
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system. The test surface of the Plexiglass v;as hand 

^ to .a flat , near-aloss . finish, 

static i^o_asur_ements a dial . indicator showed. .the 

dish to have jf/-.. 0.00 8 mm .deviation from a flat plane. 
^ho...^ssei^ly_was housed in a 1.8 m..cubical box with an 
open front. The radial flow at the disk periphery was 
ducted behind the disk by placing a 1.8 m square cover 
several boundary layer thicknesses in front of the disk 
with the test surface exposed by a large hole of slightly 
smaller diameter than the disk. 

Hot Wire System 

A single constant temperature, linearized hot wire 
was used for all measurements. Wollostan Pt-lORh 0.0025 mm 
wire with an active length of 1.0 mm was used. Calibration 
was done in the entry region of a duct flow at room temperature 
with the wire and prongs in the same orientation with 
respect to the flow as during the boundary layer measurements . 
For measurements, the wire was placed parallel to the disk 
surface with the wire axis along a radius. The wire was 
fixed in space while the disk boundary layer rotated past 
it. 

Response of a Single Hot Wire in a Three-Dimensional Disk Flow 


A single hot wire parallel to the disk surface with the 
wire axis along a radius will indicate a mean velocity having 
the magnitude of the vector composed of the axial and tangential 
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components of the flow. In non-dimensional terms, the 

ratio of this magnitude to that of the tangential component 
is given by: 




1 H 
R (G+1) 


■) 


{ 5.1) 


Provided that the term i(i j ^ij ) ^ is small, the hot wire 
will provide an excellent estimate of the actual tangential 
velocity distribution. For the current experiment, the 
minimum value of R was 125. For this case, Eq. ( 5 . 1 ) 
shows that the discrepancy between the actual tangential 
velocity and the quantity measured by the hot wire is less 
than one percent for all G>-0.95 or over approximately 
the inner 80% of the boundary layer thickness. Since the 
outer region of the boundary layer is characterized by 
extremely low velocities for which measurements are inherently 
inaccurate, no corrections were made to the measured profiles. 

To verify that a well behaved disk flow was present, the 
mean tangential velocity distribution was measured at 
several Reynolds numbers. Results are shown in Figure 2. 

The profiles for R = 251, 374 are in good agreement 

with the exact solution (see Eqs. (2.6) - (2.10)). The 
distribution at higher Reynolds numbers Reynolds numbers 
is expected to deviate from the exact solution due to the 
presence of the highly amplified stationary vortices. Even 
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then, the effect of the stationary vortices on the mean 
flow is found to be largely confined to the outer region 
of the boundary layer. 

Measurements of the fluctuating components of the flow 

growth rate of the stationary 

disk disturbances. The theoretical formulation of the 
problem assumes the same growth rate Im w for the three 
velocity components . Therefore, any arbitrary combination 
of , the components will also have the same growth rate. 

The single hot wire responds to the axial and tangential 
components of the flow. The growth exhibited by the 
hot wire can, therefore, be used for comparisons with 
theoretical predictions. The experimental growth rates 
will be compared with the theoretical predictions in 
the next Section. 

Mechanical Disk Vibrations 

Due to the large velocity gradient near the disk surface, 
any displacement of the disk surface relative to the hot 
wire will modulate the anemometer output. Low frequency 
displacements due to - the static or dynamic deviation of the 
disk from a flat plane are easily recognizable since they 
modulate the signal at a frequency corresponding to some small 
multiple of rotational frequency. Figure 3 shows that 
static deviation of the disk consists of four undulations 
in the surface. Figure 4 is a plot of the hot wire 
output at R - 457 and G = ~ 0.5 for one revolution of 
the disk. The signal was bandpass filtered in the range 
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250 < f < 600 hz , however, the envelope of the signal 
roughly corresponds to the disk surface undulations and 
thus partially accounts for the modulated output. 

Other sources of vibrations were the ball bearings, 
drive system ^and structural resonances of the disk support 
stand. These vibrations were of very small amplitude and 
were noted to occur in the same frequency range as the 
passing stationary vortices. These vibrations are critical 
when measuring disturbances in regions of low fluctuation intensity. 
To show the effect of these vibrations, the displacement of 
the center of the disk was monitored with a proximeter 
and the spectrum of these vibrations was compared to that of 
the hot wire output at R - 125 and G = -0.5. The 
vibration data was converted to apparent velocity fluctuations 
by using the calculated v(z) velocity gradient at the hot 
wire location. It was assumed that the vibrations at the 
center of the disk were very similar to those at the hot 
wire location. Results are shown in Figure 5 . (Both 
sets of data were bandpass filtered in the range 250<f<600 Hz) . 

The plots show that the hot wire output can be largely 
attributed to the disk vibrations. Since the vibrations ^ 
were in the same frequency range as the stationary vortices, 
the hot wire data at low intensities could not be used to 
indicate the amplification rate of the stationary disk 
disturbances. This also eliminated the possibility of 
experimentally determining the critical Reynolds number for 
the stationary disturbances. 
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6 . 


RESULTS AND DISCUSSION 

Critical and Transition Reynolds Niombers 
Some of the available experimental data for critical 
and transition Reynolds numbers are given in Table 1. 

It is apparent that there is considerable variation of 
the observed critical point. We believe that the 
variation can be attributed to the different measurement 
techniques used in the experiments. Using the Orr- 
Sommerfeld equation, we obtained a critical Reynolds 
number -171 which is in good agreement with the 

theoretical results of Brown [14] and Cebeci and 
Stewartson [3] but is considerably less than the 
observed values. The value of the critical Reynolds 
number for stationary vortices is significantly improved 
when the effects of Coriolis forces and streamline 
curvature are included. Our calculated critical 
Reynolds number of 287 is in excellent agreement with 
the value of 297 obtained by Kobayashi et al [12] 
using hot-wire techniques. Kobayashi et al [12] also 
performed a theoretical analysis in which some of the 
effects of Coriolis forces and streamline curvature were 
considered. They calculated a critical Reynolds number 
of 261. 

In order to correlate transition using stability 
theory, one has to know the experimental location of 
the onset of transition. The transition Reynolds 
numbers usually given for experiments (see Table 1) 
are the locations where transition is complete. 


63 



Gregory and Walker [21] showed that, for a slitted rotat- 
ing disk, the transition region is composed of two 
subregions: (i) a vortex region and (ii) an inter- 

mediate turbulent region where the intermittancy factor 
Y varies from 0 to 1. Stability theory is only 
applicable up to the point where the first turbulence 
burst appears (y = 0) . Gregory and Walker obtained 
R = 505 and 524 for y= 0 and y = 1, respectively. 

Chin and Litt [23], using an electrochemical technique, 
observed that the transition was complete at R = 592. 

They also observed that the vortices start breaking 

down into turbulence at R = 510. We believe that this 

T 

result should be taken as the relevant location for the 
onset of transition for the purposes of comparison 
with stability theory. Further evidence that the initiation 
of transition occurs at R^ 510 is provided by Kobayashi 
et al who observe that the disturbances are non-linear 
at R = 500. Usually the non-linear region is narrow 
so the onset of transition soon follows. Further, Federov 
et al [13] observed turbulent flow at R - 515. On 
the basis of all this evidence we take R^ 510 as the 
location of the onset of transition. 

Growth of Infinitesimal Disturbances 

Disturbances of all frequencies may be present in 
natural transition. We follow the evolution of several 
different modes and the one which gives the highest inte- 
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grated growth factor is used to correlate transition. 
Stationary disturbances were found to give the highest 
N factor for rotating disk flow over all positive real 
frequencies. Disturbances with negative phase velocities 
can give slightly higher N factors but they are of 
no consequence in the process of transition. 

It was shown in [6] that envelope method is a 
reliable tool for transition prediction in three 
dimensional flows. First, we report calculations using 
the Orr-Sbmmerfeld equation. The resulting N factors 
are compared with those of Cebeci and Stewartson [3] 
in Figure 6. It is evident that their method predicts 
N 2i20 at transition (R = 510) while the present 
(envelope) method gives N -22 at transition. 

Cebeci and Stewartson [3] used spatial stability 
theory and in order to remove arbitrariness among the 
parameters of equation (4.1), they imposed the condition 

8a . 

^ = 0 ( 6 . 1 ) 


65 



where =■ Im(a) / = Re(B). In order to simplify 

their computations, Cebeci and Stewartson assumed that the 

maximum growth rate at any R > R is independent of 

c 

the growth direction. This condition is not 
realistic and we believe that had their growth rates 
been maximized over all possible growth directions their 
N factor at transition would be in better agreement with 
the present predictions using the Orr-Sommerf eld 
equation. 

In Fig. 7 we plot calculated temporal growth rates 
(Im((jL>)) for stationary vortices. It can be seen that 
the inclusion of streamline curvature and Coriolis 
forces have a significant stabilizing effect. Calculations 
with only Coriolis terms (as done by Lilly [15] for 
Ekman flow) were also made. These results indicate 
that streamline curvature effects must also be included 
in order to model properly the physical problem. 

Since the instability is spatial in nature, we 
transform temporal growth rates to spatial growth rates 
a using the group velocity transformation 


Im(a3) 

Re (o) ) 
a 


( 6 . 2 ) 
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The spatiaX girowth !ra.t 0 s air© plottsd in Figuir© 8 
together with experimental data. Except for the data 
at low fluctuation intensity and for data in the 
turbulent breakdown region, the agreement with the present 
stability analysis is good. 

Integrated growth rate (N factor) results are 
presented in Figure 9. The present stability theory 
gives N ^10.6 which is close to the value N = 9 
for two dimensional flows and is in the range of values 
found for swept wings [6]. It is apparent that there 
is a very significant effect on the predicted transition 
N factor when the effects of Coriolis forces and stream- 
line curvature are included. The resulting N factors 
are much more reasonable than those obtained by conventional 
stability theory where only the Orr-Soramerf eld equation 
is solved. 

Also presented in Figure 9 are experimental results 
for the N factor. The experimental amplification rate 
of the stationary vortices was determined from the rms spectra 
of the bandpass filtered hot wire output. The rms 
voltage in a narrow band centered on the frequency (438 + 8 Hz) 
of the stationary vortices being sv/ept past the probe was calculated 
for each Reynolds number. The N factor was based on this 
rms level relative to the local disk velocity. Due to the 
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problem with disk vibrations indicated earlier, it was 
assumed that = 1 and the resultant data were shifted 
at constant Reynolds number to match the theoretical growth 
curve. The data are seen to be in a fair agreement with 
the present theory over the range 400<R<500. The significant 
deviation of the data for R <400 is attributed to disk 
vibrations. The falling off of the data for R > 500 is due 
to the highly non-linear nature of the flow in this region 
and breakdown to turbulence. 

Although shifting the level of the data because of uncer- 
tainty in the value of A^ can be questioned, the fact that the 
slope of the experimental curve (a = dN/dR) matches the present 
theory (see Figure 8) in the range 400<R<500 is very encouraging. 
With recent advances in rotating equipment technology, it may be 
feasible to build a disk drive system with low enough vibration 
amplitudes in the frequency range of interest to allow an experi- 
mental estimate of A^ to be made. Work investigating this 
possibility is now underway. 

The absolute value of the fluctuation intensity is plotted 
in Figure 10. The intensities for the range in which the growth 
rate agrees with the theory is seen to be from 0.1% to 10%. ^ 



Orientation and Number of Vortices 


In the envelope method we maximize growth rates 
Im(o)) for stationary vortices over all possible wave 
angles. We find that the vortex spirals make an 
angle of 11 . 2 ® with the negative of the direction 
of disk rotation. This is in excellent agreement with 
the experimental value of 11-14® [10-13]. 

It can be shown that the number of vortices 
is given by 


= 6 R (6.3) 

where 3 is defined in (2.19) . 

Gregory et al [10] observed about 30 vortices 
in the range of Reynolds numbers range 430-530. At 
R * 430/ we obtain n = 0.0698 x 430 ^30, which is in 
excellent agreement with the experimental observation. 

If the number of vortices is to remain constant then 
3 should vary as 1/R, Our calculations do not show 
this behavior. Instead, 3 remains almost constant so 
that n varies with R. 

The continuous variation of n can not be justified 
physically. However, it is possible that the number of 
vortices does vary but not in a continuous fashion. 

There is some experimental evidence of bifurcation in which 
n undergoes discrete jumps. Fedorov et al [13] observed 
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30 vortices at R = 387. However at lower Reynolds numbers 
(R = 180 - 245) , they observe 14-16 vortices. Linear 
theory is unable to predict such bifurcation phenomena 
so nonlinear theory may be needed. 

Parallel or Type II Instability 

Lilly [15] presented numerical solutions of the 
Ekman layer problem and included the effect of Coriolis 
forces in his analysis. He found that at very low Reynolds 
numbers an instability mechanism exists whose disturbances are 
different from the stationary disturbances described 
previously. Lilly called this "parallel instability" and 
suggested that it is of viscous type since it vanishes 
at high Reynolds mmbers . He found that the critical 
Reynolds number for these fast moving disturbances is 55 
and the resulting modes are oriented at small negative angles. 
The orientation angle at the critical point is -23° 
which decreases in magnitude as the Reynolds number 
increases. A similar instability mechanism was detected 
in the experiments of Faller and Kay lor [16] (who called 
it a type II instability) and Tatro and Mollo-Christensen 
[17] . 

In our calculations, we also find travelling 
disturbances (f~ 100 Hz relative to the disk) at low 
Reynolds numbers. The critical Reynolds number for these 
disturbances is calculated to be 49. The critical 
parameters are 

ct = 0.27 3 = - 0.137 = 0.146 
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(6.4) 



This corresponds to a wave oriented at an angle of -26.9°. 
These disturbances have much lower grov/th rates than the 
stationary ones . 

CONCLUSIONS 

The growth of instabilities in the three dimensional 
flow due to a rotating disk is studied both theoretically 
and experimentally. The experiments show clearly a 
region of linear grov7th that is in good agreement with 
linear stability theory that includes the effects of 
Coriolis forces and streamline curvature. Therefore, 
the e^ method gives good results for transition 
prediction in these three dimensional boundary layers 
with N of the order 11. 
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Table 1. Critical and Transition Reynolds Numbers for Rotating Disk Flow 


Reynolds Number 

Investigators Critical Transition Unset o± Method of 

Transition Investigation 
• (estimated) 


Smith 

[9] (1947) 

460 

557 

- 

hot-wire probe 

Gregory et al 
[10] (1955) 

430 

530 

- 

visual 
(China clay 
technique) 

Cobb & 
[23] 

Saunders 

(1956) 

447 

490 . 

- 

heat transfer from 
the disk 

Gregory & Walker 
[21] (I960) 

367 

524 

505 

acoustical 

slitted disk 

Chin & 
[2:j 

Litt 

(1972) 

412 

592 

510 

mass transfer 

coefficient using 

electrochemical 

technique 

Fedorov et al 
[13] (1976) 

387 

515 

- 

visual (Napthalene) 
acoustical 

Clarkson et al 
[11] (1980) 

532-621 

562-680 

_ 

visual 

(dye in water) 

Kobayashi et al 
[12] (1980) 

297 

566 

500 (non- 
linear 

oscillations ) 

hot wire probe 

Present results 

287 

- 

- 

calculations using 


stability theory 
includina Coriolis 
force and strsamlin 
curvature effects 
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FIGURE CAPTIONS 


Figure 1. 
Figure 2. 

Figure 3. 
Figure 4, 
Figure 5. 

Figure 6 . 
Figure 7. 


Disk system lay out — top view 

Normalized mean tangential velocity profiles for the 
KArmin rotating disk flow 

Measured deviation of disk surface from a flat plane 

Time variation of fluctuation intensity 

Comparison of disk vibration spectrum with hot wire 
output. (a) Hot wire output (b) Proximeter 

Integrated growth factor using Orr-Soramerf eld equation 

(a) Present calculations 

(b) Cebeci and Stewartson [3] 

Temporal growth rates for stationary vortices 

(a) Orr-Sommerf eld equation 

(b) Orr-Sommerf eld equation with Coriolis force effects 
included 

(c) Orr— Sommer f eld equation v;ith Coriolis force and 
streamline curvature effects included 


Figure 8. Spatial growth rates for stationary vortices 

(a) Orr-Sommerf eld euqation 

(b) Present theory 
(+) Experimental data 


Figure 9. Integrated growth factor for stationary vortices 

(a) Orr-Sommerf eld equation 

(b) Present theory 
(+) Experimental data 

Figure 10. RMS intensity of velocity fluctuations as a function of 
Reynolds number 
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Comparison of disk vibration spectrum with hot wire 
output. (a) Hot wire output (b) Proximeter 


Figure 5. 
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V. 


Subcritical Transition to Turbulence in Plane Channel 
Flows 


While experiments^ show that incompressible plane 
Poiseuille and plane Couette flow may undergo transition 
to turbulence at Reynolds numbers R of order 1000, 
linear stability analysis of these plane parallel flows 

gives critical Reynolds numbers of 5772 for 

2 3 

plane Poiseuille flow and ~ for plane Couette flow . 

This discrepancy between theory and experiment suggests 

that the mechanism of transition is hot properly represented 

by parallel-flow linear stability analysis. In this Letter, 

we present a new linear three-dimensional mechanism that 

predicts transition at Reynolds numbers in good agreement 

with experiment for both plane Poiseuille and plane Couette 

flows. Here we present the theory applied to plane 

Poiseuille flow, defined as flow between fixed parallel 

plates that is driven by a pressure gradient. 

We begin by studying two-dimensional travelling-wave 

solutions to the Navier-Stokes equations : 

v(x,z,t) = F(x-ct,z) 

where c is a real wave speed, x is the downstream 

coordinate and z is the coordinate perpendicular to the 

channel walls at z = ±1. No-slip boundary conditions are 

applied at the walls and 2Tr/a periodicity in x is assumed. 

2 

For all R , one solution is the laminar flow (1-z ) , where 


( 1 ) 
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R = 1/v and v is the kinematic viscosity. For R >2900, 

up to 2 other solutions (neglecting an arbitrary phase 1 
may exist for any given a The locus of points in 
(E,R,a) space for which these solutions exist is called 
the neutral surface. Here E is the energy of the flow 
relative to that of the l^inar flow. A slice of the 
neutral surface for given subcritical Reynolds number 
(290C< R < 5772) is shown in Fig, 1,^ 

If a one-dimensional phase space representation were 
appropriate to describe the behavior* of flows off the neutral 
surface, E would evolve according to 

( 2 ) 

Typically the critical points of (2) are alternately stable 
and unstable, so the lower branch (LB) solutions on the 
subcritical neutral surface plotted in Fig. 1 are unstable 
while the upper branch (UB) solutions are stable. 

While these stability predictions are correct, the 

evolution of two-dimensional flows is not restricted to a 

one-dimensional phase space. Projections of numerical 
5 

solutions of the two-dimensional Navier-Stokes equations 
on the two-dimensional phase space (/e^, are plotted 

in Fig. 2. Here E^ is the kinetic energy in that part 
of the flow that depends on x like Orbits of 

solutions with initially large energies do not follow simple 


dE ^ 
dt 


f CE) 



curves. The time-dependent evolution of two-dimensional 
flows evidently requires a multi- (likely infinite) dimensional 

g 

phase space. Thus, Landau-Stuart-Watson nonlinear 
stability theory, which gives evolution equations of 

- - . . - . - - 7 

the form (2) can not be valid away from the neutral surface. 

Several other features of Fig. 2 are noteworthy. First, 

the two orbits in the lower left hand corner illustrate the 

existence of a threshold energy (near that of the LB 

solution) below which disturbances decay. Second, solutions 

with energies less than that of the UB solution (indicated 

by the point marked 'steady solution* in Fig. 2) can overshoot 

the UB energy by factors of 4 or more. Third, and most 

importantly, typical solutions quickly evolve to a state 

within a band of quasi-equilibria and, then, only very 

slowly approach the steady UB solution. The time scale 

for initial adjustment to a quasi-equilibrium state is 

of order the eddy circulation time 1//E, (i.e., of 

order 10) while the time scale for approach to the equilibrium 

state is of order the diffusion time 1/v (i.e., of 

order 1000 - 10000) . In the quasi-equilibria, the spanwise 

8 

vorticity must be nearly constant on streamlines , so that 
equilibrium is achieved by diffusion of vorticity. In fact, 
vorticity can vary by at most 0(v) along interior stream- 
lines of the equilibrium flows. Nearby flows must have the 
same property implying the existence of quasi-equilibria 
evolving only on, a diffusive time scale. 

The quasi-equilibria are the basis of our transition 

mechanism in plane Poiseuille flow as direct numerical 

9 

solution of the Navier-Stokes equations shows that they 
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are strongly unstable to infinitesimal three-dimensional 
disturbances. In Fig, 3, we plot the evolution of (initially 
small) three-dimensional disturbances superposed on finite- 
amplitude two-dimensional motions. Evidently, the three- 
dimensional disturbances quickly achieve a form that grows 
exponentially in time for R>1000, The growth rate of the 
three-dimensional disturbances is rapid with their amplitude 
increasing by a factor of about 10 in a time of 10. 

This short time scale for subcriticai three-dimensional 

growth should be contrasted with the long time scale of 

~ 2 

order 1000 for evolution of supercritical Orr-Sommerf eld modes. 

There is strong evidence that this instability is a 
physically relevant one in that it is fairly insensitive to 
initial conditions and has small threshold energies. It is 
necessary to distinguish here between this instability and 
the ensuing transition to turbulence. If the two-dimensional 
flow persists sufficiently long for the three-dimensional 
perturbations to attain a finite amplitude, direct numerical 
simulation has shown^ that the resulting three-dimensional 
flow quickly develops a turbulent character with strongly 
non-periodic behavior. Thus to 'predict’ transition one 
must know the initial two-dimensional and three-dimensional 
energies as well as their respective time scales. For instance, 
the most dangerous three dimensional instability for given 
two-dimensional energy is not necessarily the most likely to 
force transition if the two-dimensional state is outside the 
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band (in wavenumber) of quasi-equilibria. It is possible to 
use our methods to construct a neutral surface for transition 
in any given (presumably large) parameter space. However/ 
we confine attention here to demonstrating that our mechanism 
predicts transitional Reynolds numbers in accordance with 
experiment. 

The exponential growth illustrated in Fig, 3 suggests 

10 

that a linear instability mechanism is involved. Assuming 
a flow of the form 

vCx,t) = F(x-ct,z) + e Re[GCx-ct^z)e^^‘*'^^^] , (3) 

substituting into the Navier— Stokes equations, and linearizing 

with respect to e , a linear eigenvalue problem for q 

results. The Galilean transformation to a reference frame 

moving with the phase speed c eliminates time-dependent 

coefficients, so the problem is separable in t. The resulting 

eigenvalue problem has been solved numerically using Chebyshev 

polynomial expansions in z and highly truncated Fourier 

series expansions in x for F and G. In Fig. 4, we plot 

11 > 

the maximum growth rate , ReCcJ) f vs the spanwise wavenumber 
3 for R = 4000, a =1,25, The results of direct numerical 
simulations fcf. Fig, 3) are also plotted in Fig, 4. Evidently, 
the large growth rates observed in the direct numerical simulations 
can be explained by this linear eigenvalue problem. 

Note that the linear theory presented above can be 
extended to Reynolds numbers below 2900 by freezing the 
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quasi-equilibria which evolve very slowly compared to the 
rapid exponential growth of the three-dimensional perturbations. 

For R>1000, the quasi-equilibria decay sufficiently slowly 
that three-dimensional perturbations can grow, overwhelm the 
two-dimensional flow, and break down to turbulence. 

The rapid growth rates described above are due to the 

combined action of vortex stretching by the nearly inviscid 

two-dimensional steady motion F and tilting of the vortex 

lines of F by the perturbation G. By itself, vortex 

stretching by F can not give exponential growth rates 

1 0 

because of the tv70-dimensional anti-dynamo theorem. 

Detailed flow visualizations of the instabilities described 
here will be given elsewhere. It will be shown that 
three-dimensional perturbations grow on a time scale of 
order I/i/e^.q/ which must be shorter than the decay 
time of the two dimensional motion for the instability 
to be effective. The sharp cutoff in growth rate Re (a) for 
small 3 observed in Fig, 4 reflects a threshold of streamwise 
vorticity for stretching to persist. 

. 9 

Direct numerical simulations of transition in plane 
Couette flow show that while there is no evidence that 
equilibria of the form (1) exist, the three-dimensional ^ 

instability process outlined above is still effective down to 
Reynolds numbers of order 1000. While there are no quasi- 
equilibria in plane Couette flow that evolve on purely 
diffusive time scales, the decay rates of finite-amplitude 
two-dimensional disturbances are still several eddy circulation times 
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This implies that the threshold three-dimensional 
energies in plane Couette flow are somewhat larger than in 
plane Poiseuille flow. However, the resulting instability 
is at least as strong and turbulence quickly ensues. 
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Figure Captions 

Fig. 1. 

A subcritical (E,a) slice of the neutral surface for 
plane Poiseuille flow at R = 4000, The stability of 
solutions is indicated by the arrows. The behavior shown 
in this plot is typical for 2900 <R< 5772. 

Fig. 2 . 

A phase portrait of disturbances to laminar plane 
Poiseuille flow in space at R = 4000, a = 1.25. 

The dots, equally spaced by 1.25 in time, indicate the 
evolution of perturbations from initial conditions proportional 
to the least stable Orr-Sommerfeld eigenfunction at this 
(a,R). Note the existence of a band of quasi-equilibria. 

Fig. 3. 

A plot of the growth of three-dimensional perturbations 
on finite-amplitude two-dimensional states in plane 
Poiseuille flow at (a, 3) = (1.32, 1.32). Here ^2-D 
the total energy (relative to the laminar flow) in wave- 
numbers of the form (na,0) , while ^ 3 -p the total 

energy in wavenumbers (na,3) • For R > 1000 we obtain 
growth and for R = 500 decay. The growth rate of the 
three-dimensional disturbance amplitude at R = 4000 
is about 0.18 and depends only weakly on R 

for larger R. The initial conditions are superpositions 
of the laminar flow, a [large (E 2 _j^— 0 . 04 ) ] two-dimensional 
Orr-Sommerfeld mode with wavevector (ot,0) and a [very 
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small ^^3— D ~ ^ ^ ■thr©©’^dimensional traricvoa^a© 

Orr-Sommerfeld mode with wavevectors (0,3). 

Fig, 4. 

A plot of the growth rate a of three-dimensional 
perturbations as a function of 3 at R = 4000, a = 1.25, 
Note the good agreement between the linear calculation and 
the 2-mode direct simulations. Increasing the number of 
retained modes in x increases the growth rates. However, 
the error in the 2-mode model is not large. 
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Fia 3 A plot of the growth of three dimensional 

perturbations on finite-amplitude two-dxmensronal 

States in plane Poiseuille flow at (a , g) - (1 . 32 , 1. . 

Here E^n i" the total energy (relati.^ to the lamrnar 
nL) in'Savenuirbers of the form (na,0), whrle 13 .^ 
is the total energy in wavenumbers 
R>1000 we obtain growth and for R ~ decay. 

■ growth rate of the three-dimensional dist^ance 
amplitude at R = 4000 is about 0.18 ^ 

depends on ly weakly on R for larger R. The 


3-D 
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initial conditions are super- 
positions of the laminar flow, 
a [large (E2 -d= 0»04)1 two- 

dimensional Orr— Sommerf eld 
mode with wavevector 
(a,0) and k[very small 
(E3 _j 3=10“16) ] three- 
dimensional transverse 
Orr-Sommerf eld mode 
with wavevector (0, 3 ) . 


Fig, 4. A plot of the growth rate a of 
three-dimensional perturbations as a function 
of 3 at R= 4000, a = 1.25. Note the good 
agreement between the linear calculation and 
the 2-mode direct simulations. Increasing the 
number of retained modes in x increases the 
growth rates. However, the error in the 2 -mode 
model is not large. 



VI . Finite Amplitude Stability of Axisyimnetric Pipe Flow 


abstract 


The stability of pipe flow to axisyinmetric disturbances 
is studied by direct nximerical simulation of the incompressible 
Navier-Stokes equations. There is no evidence of finite- 
amplitude equilibria at any of the v/avenumber/ Reynolds 
number combinations investigated, with all perturbations decaying 
on a time scale much shorter than the diffusive (viscous) 
time scale. In particular, decay is obtained where 
amplitude -exp an si on perturbation techniques predict 
equilibria, indicating that these methods are not valid 
away from the neutral curve of linear stability theory. 



1. INTRODUCTION 

It is generally agreed that pipe Poiseuille flow 
(also called Hagen-Poiseuille flow) is linearly stable to 
all disturbances (both axi symmetric and non- axi symmetric) 
at all Reynolds numbers (Sexl 1927, Lessen et al 1968, 

Davey & Drazin 1969, Metcalfe^ Orszag 1973, Salwen et 
al 1980) . Therefore, the explanation of the observed 
transition to turbulence in this flow requires finite- 
amplitude instabilities. 

Finite-amplitude stability analyses of pipe flow 
have so far been restricted to axi symmetric disturbances 
(Davey & Nguyen 1971) and even these results are not 
without controversy (Itoh 1977, Davey 1977). In this 
paper, we use spectral methods to investigate numerically 
the behavior of finite-amplitude axisymmetric disturbances in 
pressure-driven pipe flow. The basic question concerns the 
existence of finite-amplitude equilibrium states of this 
flow. If such states do not exist then pipe flow is stable 
to all axisymmetric distrubances . Available finite- 
amplitude analyses predict equilibria, however they are 
in disagreement over both results and methodology. 
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2 . NUMERICAL METHODS 


The axi. symmetric incompressible Navier-Stokes 
equations are, in rotation form. 


9u 


4 12 

+ VOJ + i (D^ + U 


9x R 


( 1 ) 


9v 

3t 


- U03 + I (d2 - ^ + 3 )V 

r 9x 


( 2 ) 


^Uu) + (rv) = 0 


(3) 


where u and v are the velocities in the x and r 

. 2 1 
directions, respectively, D = — 9/9r(r9/9r), and 

0 ) = 3v/9x “9u/3r is the azimuthal vorticity. 

The boundary conditions on the velocities are 


u, v/r bounded (r = 0) , u,v = 0 (r = 1) (4) 

R is the Reynolds number based on pipe radius 
and center-line velocity. The constant pressure gradient 
term is assumed to be that of the laminar flow, namely 
4/R, and R is the disturbance pressure head. 

We discuss briefly below four features of our numerical 
methods: spectral representation; time-stepping; 

operator inversion; and initial conditions. The major 
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difference between the present pipe flow calculations 
and our plane channel flow simulations (Orszag & Kells 
1980; Patera & Orszag 1980) is that variable-coefficient 
equations must now be solved implicitly, whereas 
in planar geometry only constant coefficient equations 
require’ implicit solution. 

The velocities are expanded as 

u(x,r,t) = I 

ln| <N 

v(x,r,t) = I 

|n|<N 


p=0 


u(n,p,t)e ^2p^^^ 


(5) 


I v(n,p,t)e^°^^^T^ _ , (r) (6) 


p=0 


2p+l 


where the Chebyshev polynomials defined by 


T (cos 0 ) = cos q 0 


The even (odd) nature of u(v) follows naturally from the 
axisymmetry of the problem. Boundedness at the origin 
is ‘then automatically imposed. Periodic boundarv. . 
conditions are applied in x with periodicity interval 
2ir/a. 
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A fractional time stepping method (Orszag & Kells 
1980) is used, each full step consisting of (i) 
a non-linear step, (ii) a pressure step, and (iii) 
a viscous dissipation step. For the first fractional 
step a second order Adams-Bashforth method is used: 


/vn+1 n , . . ,3 n n 1 n-1 n-1 , . 

u =u+At(jV0) “'2'^ ^ ^ 4 /R) 


''n+1 n , , 3 n n , 1 n-1,^ n-1, 

V = V + At (-• 2^ GO + ^u 0) ) 


where the superscript refers to time step. The non- 
linear terms are calculated efficiently using transform 
methods and collocation (pseudospectral) techniques 
(Gottlieb & Orszag 1977) . Products are evaluated in 
physical space while derivatives are calculated in 
spectral space. Transformations between the two 
representations are done using extensions of the 
fast Fourier transform algorithm. In the remainder of 
this Section it 1^=^ assumed" that the vslocities are in 
mixed representation, Fourier in x. but physical in 
r. The axial wavenumber of a Fourier mode will be 
denoted by y . 
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Next, incompressibility is imposed with the 
pressure step 


^ n+1 ^ -n+1 


- n+1 

V - 

itru^^'^^+ 


/\ 

V 


n+1 


_9_ 

3r 


- iyAtn* (r < 1) 

* 

- At (r < 1) 

(rv^"^^) =0 (r :< 1) 


(7) 

(8) 

(9) 


$ n+1 


0 (r - 1) , 


( 10 ) 


* 

where we write n ,5 rather than II to indicate that the 
pressure obtained here is an intermediate result. Eqns. {7)>-(10) 
can be combined to give a single equation for H , 


(D^ - 


Y^)1T 


= XY u 


n+1 


+ i 

r 


9r 


(rv"+l) 


Cr < 1) 


Cll) 


an 
a r 


= 0 (r = 1) . 


( 12 ) 


Note that Yi^ is e::panded in an even series of Chebyshev 
polynomials like (5). The solution of (11) -(12) is 
discussed below. 
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In the final fractional step, viscous effects are 


included using the Euler backward scheme 


,n+l ?n+l , At 2 . n+1 , 

= ^ + — (D “Y )u (r< 1) 


u 


n+1 


= 0 (r = 1) 


n+1 ^n+1 At 

r = V 


(r<l) 


= 0 Ir = 1) . 


(13) 

(14) 

(15) 

(16) 


The overall scheme is only first order accurate in time 
because the viscous and pressure operators do not commute 
Higher order accuracy in time may be obtained by extra- 
polation methods . 

The implicit parts of the procedure given above all 
involve solution of an equation of the form 


(L - 3 ^) (j) = f (r < 1) 


(17) 


a({> + 



r = 1 


(18) 
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for each Fourier mode, where L is a second order 

2 

(Laplacxan) operator in r, 3 depends only on the 
x-Fourier index (not r) , and a and b are constants 
independent of Fourier mode. We discuss briefly here 
the discretization of L and the method of inversion. 

For planar geometries L can be written as a tri- 
diagonal system using a Chebyshev tau -method (Orszag 
and Kells 1980) . This system can then be inverted in 
0(P) operations for each x-Fourier mode. The 
curvature terras introduced by the cylindrical geometry 
destroy the tridiagonal property of the tau-method 
matrix, and collocation then becomes more attractive due 
to the ease with which variable coefficient problems 
can be handled. 

To solve the full matrix equations resulting from 
the collocation approximation of (17) - (18) , an 
eigenfunction solver is used that reduces the operation 

3 2. 

count from 0 (P ) to 0 (P ) while only requiring the 
storage of one P x p matrix for given (L,a,b) . More 
precisely, we diagonalize the collocation approximation to 
L as 

L = A'i' 

The solution to (17) - (18) can then be written as 
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<ji = 1'“^ (A - T f 


The diagonalization (independent of Fourier mode) is 
done in a pre-processing stage. 

Finally, the initial conditions for the runs 
presented here are of the form 


v(x,r,t = 0) = 


2 

(l-r )x + A v_ (x,r) 

Li 


where v is an eigenfunction of the fourth-order 

1j 

linear stability equation obtained from (1) - (3) by 
assuming a solution of the form 


V = (l-r2)3 + 0^(r) 


and linearizing with respect to e. The magnitude of the 
Perturbation is characterized by its energy relative to 
that of the unperturbed flow ;• 


E - 12 / (u + V ) r . dr 

0 


C19) 


The details of the linear problem are well -documented 

(Lessen et al 1968, Davey & Drazin 1969, Salwen et al 
1980) and will not be discussed here. The numerical 
procedure used to determine the eigenvalue co and 
eigenfunction i|;(r) for given R, a is similar to that 
described by Orszag (1971) for planar geometries, except 
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that collocation rather than the tau-method is used. 


3 . RESULTS 

Before investigating finite-amplitude behavior it 
is necessary to confirm that the direct simulation 
gives decay rates and phase velocities in agreement 
with those predicted by linear theory. The results 
of two such tests are summarized in Table 1, where 
it is seen that the code adequately resolves both center 
and wall modes at modest time steps^ (even without 
using extrapolation methods to reduce the first order 
error in time) . 

The results of linear theory can also be used to 
test whether interactions between a primary mode and 
its harmonic are accurately simulated. For pipe 
flow, a center mode with wavenumber ot nearly resonates 
with its haimionic in the sense that the phase speed of 
the mode with wavenumber 2a is very close to that 
of the primary. If one assumes they resonate exactly 
(i.e. , the harmonic will obey a forced 

amplitude equation of the form 



and thus grows secularly in t for short times and 

attains a maximum at 



iln(-2(o^^) -iln(-a32^) 


-2o)ii + 


( 20 ) 


As the modes are not exactly phase-locked, we would expect 

* * 

the actual maximum to occur at t < t . This behavior 

r 

* * 

is verified in Fig. 1 by a plot of t /t^ at a = 1 
for various Reynolds numbers. The maximum deviation 
between (20) and our direct simuJ ations, is -3%. 

Next, we study finite-amplitude disturbances 
predicted to be dangerous by the method of false problems. 
Davey & Nguyen (1971) find that the two disturbances 
tested above at very small amplitude (see Table 1) have 
threshold (unstable equilibrium) energies (19) of 
E - 0.003. The method of Itoh (1977) as applied by 
Davey (1978) indicates that the center mode should 
decay at finite amplitude , however it too predicts a 
small- amplitude equilibrium state for the wall mode. 

To test these theoretical results, the .same series 
of runs reported in Table 1 were repeated except that the 
axial and radial resolution was increased, the time 
integration was taken to a final tim.e of T = 20 rather 
than T = 10, and the initial energies of the disturbance 
were taken to be 0.01. The results for the wall mode 
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and center mode are shown in Figs, 2 and 3 respectively, 
as plots of the logarithm of the primary and secondary 
energies as a function of time. There is apparently 
no evidence of equilibria. Runs at lower and higher 
initial energies (e.g, E - 0,04) decay in a similar 
manner. 

The lack of equilibria reported above does not 
preclude their existence for other Reynolds number/ 
wavenumber combinations. However, in a variety of runs, 
we have found no finite-amplitude steady-states. The 
results of a typical run at R = 4000 , a = 1,0, co = 0.3783 - 
i 0.1025, are plotted in Fig, 4. From Fig. 4 we 
infer that the disturbance at R = 4000 decays in a 
time very short compared to a diffusive time scale, and 
is therefore consistent with the absence of equilibria 
(Orszag & Patera 1980) . 

The time scale for decay of finite-amplitude axi- 
s^nranetric states is central to an understanding of three- 
dimensional transition. The three-dimensional instability 
mechanism leading to transition in plane Poiseuille and 
Couette flows (Orszag & Patera 1980) develops on a time 
scale of order 1 //e, where E is the energy of the two- 
dimensional disturbance. For the disturbance plotted in Fig 
1 //e is significantly shorter than ' the ‘tXme scare“on~whioh^^ 
the perturbation decays. 
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Thus, it seems that the instability we found in planar 
channel flov/s may be relevant to transition in pipe 
flow as wel], . This possibility will be investigated 
further in a later paper in which the behavior of 
non-axisymmetric disturbances to decaying axisymmetric 
states will be considered. 

Our results indicate that the method of false 
problems is not a valid procedure for investigating 
finite-amplitude axisymmetric perturbations to pipe 
flow. We do not attempt a critique of these methods 
here except to emphasize a point made by Herbert (1977). 
Herbert commented that the retention of only the first 
term in the amplitude expansion of the frequency without 
knowing the convergence properties of the series can lead 
to incorrect conclusions, especially in cases (such as 
pipe flow and plane Couette flow) where there is no linear 
neutral curve. The radius of convergence of the 
amplitude expansion may simply be too small to predict 
equilibria. For example, numerical simulations of plane 
Couette flow (Orszag & Kells 1980, Patera & Orszag 1980) 
do not confirm the existence, of two-dimensional equilibria 
predicted by Davey & Nguyen (1971) on the basis of amplitude 
expansions. The direct iteration procedure of Herbert 
(1977) ' bypasses this problem and can predict the existence 
of equilibria as well as their threshold energy. 
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However, a limited survey of the available phase space 
has not yet yielded any f inite-ampli tude solutions. 

We suspect there are none. 
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Table 1. Behavior of Linear Modes 



wall Mode 

Center Mode 

R 

1600. 

500. 

a 

00 

to 

6.2 

Re w 

1.5849 

5.8852 

Im ^ 

-0.5396 

-0.3917 

Perturbation Energy 
Spatial Resolution 
(2NX (P+1) ) 

1x10-10 

ixio’^^ 

8x33 

8x17 

At 

0.005 

0.01 

Final Time , T 

10. 

10. 

Computed Re^ w 

1.5836 

5.9146 

Computed Iin 03 

-0.5410^ 

-a.3876 
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FIGURE CAPTIONS 


Fig. 1 The ratio of the tine at which the harmonic 
attains its maximum to the time predicted by linear theory 
(assuming perfect phase-locking) is plotted as a function 
of Reynolds number when a = 1. The actual (computed) 
time is less than the predicted linear theory time 
because the real frequency of the harmonic is not 
exactly twice that of the primary. 

Fig. 2 Decay of a wall mode at R = 1600, a = 5.8 from 
an initial energy ( E = 0.01), larger than the equilibrium 
value predicted by the method of false problems. Higher 
energy disturbances also decay. Hera E is the energy 
of the disturbance relative to the mean flow [defined 
in (19)]. Here N = 8 and P = 64 in (4)- (5). 

The accuracy of this and other runs was tested by changing 
N, P and the tim.e step At. 

Fig. 3 Decay of a center mode at R = 500, a — 6.2 from 
an initial energy (E = 0.01) larger than the equilibrium 
value predicted by perturbation theory. Here N = 8 and 
P = 32 in (4)- (5) . 



119 



Fig. 4 Decay of a wall mode at R=4000, a=l. The 
decay occurs on a time scale much shorter than the 
diffusive scale indicating the absence of equilibria. 
Here N = 8 and P==32 in (4) - (5). 
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1.5 


Fig. 1 The ratio of the tine at which the harmonic 
attains its maximum to the time predicted by linear theory 
(assuming perfect phase-locking) is plotted as a function 
of Reynolds number when a = 1. The actual (computed) 
time is less than the predicted linear theory time 


because the real frequency of the hax-monic is not 


exactly twice that of the primary. 
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3 Decay of a center mode at k - 


„ (F = 0 01) larger than the equilibrium 
initial energy (E u.ux; 

T.T — R and. 
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Fig. 2 Decay of a wall mode at R - 1600, a = 5.8 from 
an initial energy ( E = 0.01), larger than the equilibrium 
value predicted by the method of false problems. Higher 
energy disturbances also decay. Here E is the energy 

V 

of the disturbance relative to the mean flow [defined 
in (19)]. Here N = 8 and P = 64 in (4)- (5). 

The accuracy of this and other runs v/as tested by changing 
N, P and the time step At. 



j R — 500 , a = 6.2 from 

3 Decay of a center mode at R - 50U, 

inltl.l e„.r„ (= - «-"ll 

ue predicted by perturbation theory. Here N 8 a 


P = 



124 




1 . Report No. 2. Government Accession No. 

NASA CR-165661 

3. Recipient's Catalog No. 

4. Title and Subtitle 

ADVANCED STABILITY ANALYSES FOR LAMINAR 
FLOW CONTROL 

5. Report Date 

February 1981 

6. Performing Organization Code 

7. Author(s) 

Steven A. Orszag 

8. Performing Organization Report No. 

CHI Report No. 43 

10. Work Unit No. 

9. Performing Organization Name and Address 

Cambridge Hydrodynamics, Inc. 
P.O. Box 249, M.I.T. Station 
Cambridge, Massachusetts 02139 

11. Contract or Grant No. 

NASl-15894 

13. Type of Report and Period Covered 

Final 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20545 

14. Sponsoring Agency Code 


15. Supplementary Notes 


Langley Technical Monitor: Dennis M, Bushnell 


16. Abstract 


Five classes of problems are addressed, (i) The extension of the SALLY 
stability analysis code to the full eighth-order compressible stability 
equations for three-dimensional boundary layer flows is discussed. The 
resulting computer code seems to be at least an order of magnitude 
faster than previous codes, (ii) A comparison of methods for prediction 
of transition using SALLY is given for incompressible flows. It is 
shown that the envelope method is most efficient, (iii) A study of 
instability and transition in rotating disk flows is given in which the 
effects of Coriolis forces and streamline curvature are included. Good 
agreement with experiments performed at NASA LaRC is demonstrated. 

(iv) A new 'linear* three-dimensional instability mechanism that 
predicts Reynolds numbers for transition to turbulence in planar shear 
flows in good agreement with experiment is introduced and studied. 

(v) A study of the stability of finite-amplitude disturbances in 
axisymmetric pipe flow shows the stability of this flow to all nonlineai 
axisyimnetric disturbances, in contrast to the predictions of nonlinear 
stability theory. 


17. Key Words (Suggested by Author(s)) 

Stability theory Numerical methoc 

Laminar flow control 
Nonlinear stability 
Transition Boundary layers 

Shear flow Spectral methods 

18. Distribution Statement 

is 

Unlimited - Unclassified 

19. Security Classif. (of this report) | 

Unclassified | 

20. Security Classif. (of this page) j 

Unclassified 

21. No. of Pages 

124 

22. Price* 


N-305 


For sale by the National Technical Information Service, Springfield, Virginia 22161 




















